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ABSTRACT 


The inverse problem of estimating the open channel flow roughness (Man- 
nings n) is solved using an optimization algorithm. Measurements of flow variables 
(depth of flow and discharge) at different locations and times are used as inputs to the 
optimization model. Unsteady flow conditions are considered. Both, single value of 
Mannings “n” throughout the channel length and multiple values of “n” at different 
channel reaches are considered. The Sequential Quadratic Programming algorithm 
is used to solve the optimization model. The optimization model embeds the finite 
difference forms of the governing equations of open channel flow. The Preissmann 
scheme is used for discretization of the flow equations. 

The performance of the developed optimization model is evaluated for dif- 
ferent scenarios of data availability. Solution results of the model for illustrative 
problems establish the potential applicability of the proposed approach. 
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Chapter 1 

INTRODUCTION 


1.1 Introduction 

An accurate estimation of the roughness coefficient is of prime importance to any study 
involving open-channel flow such as: 

• Flood estimation, routing and damage mitigation. 

• Optimal design and operation of canal and surface irrigation systems. 

• Estimation of degradation and aggradation of channels due to natural and human- 
made causes and 

• Other river engineering studies. 

Mathematical models which are routinely used in these studies require calibra- 
tion with regards to the roughness coefficient before they are applied. In fact, the 
whole process of effective utilization of surface water resources requires a knowledge of 
the roughness coefficient. 

Estimation of the roughness coefficient or the Mannings “n” for a natural 
channel is not a trivial task because it depends on various factors such as: 
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1. surface roughness. 

2. vegetation. 

3. channel irregularity. 

4. obstructions. 

5. channel alignment. 

6. sedimentation and scouring. 

Many empirical procedures (French, 1990) have been suggested in the past for 
estimating the value of n. The Soil Conservation Service (SCS) method (Urquhart 
1975) for estimating “n” involves the selection of a basic “n” value for a uniform, 
straight and regular channel in a native material and then modifying this value by 
adding correction factors. These correction factors are determined by a critical consid- 
eration of some of the factors mentioned above. The selection of correction factors in 
this method is highly heuristic. 

A second method of estimating “n” for a channel involves the use of tables 
of values. Chow (1959) presented an extensive table of “n” values for various types 
of channels. In this table, a minimum, a normal and a maximum value of “n” are 
stated for each type of channel. This method gives satisfactory results for human- 
made channels. However, it is not suitable for estimating the roughness coefficient of 
the natural channels. 

It is the belief of the U.S. Geological Survey that photographs of channels of 
known resistance together with a summary of the geometric and hydraulic parameters 
which define the channel for a specified flow rate, can be useful in estimating the 
resistance coefficient (Barnes 1967). The U.S. Geological Survey maintains a program 
which trains engineers in the estimation of channel resistance coefficient. It is claimed 
that the trained engineers can estimate resistance coefficeints with an accuracy of ±15% 
under most conditions. However, it should be noted that this is not a formal procedure 
and a satisfactory application is highly dependent on the person’s previous experience. 
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From a theoritical viewpoint, a value of the resistance coefficient can be es- 
tinaated from the velocity and the depth measurements because these depend on the 
roughness. Chow (1959) suggested the following equation for estimating “n” from the 
velocity measurements at a cross-section. . 


(rr — l)/i6 
6.78(a: -h 0.95) 


( 1 . 1 ) 


where: 


Uo. 2 =Velocity at of the depth 

Uo.8=Velocity at of the depth 

h=depth of flow 


X = 


Uo.s 


( 1 . 2 ) 


The above equations (1.1 and 1.2) are applicable only when channel is wide and 
the velocity measurements are made for uniform flow conditions. However, the idea 
of estimating the Mannings “n” from the field measurements of discharge and depth 
i.e. the “inverse problem” can be used for non-uniform flow conditions in irregular 
channels. The Standard Step Method (Chow 1959, Subramanya 1984) can be used 
to simulate the steady, non-uniform open channel flow. The Mannings “n” for the 
channel may be determined by a trial and error procedure such that the difference 
between the simulated and the measured water surface profile is minimum. Two points 
should be noted here with regards to the inverse problem. Firstly, a formal search 
technique (Rao, 1984) may be used for determining the Mannings “n” instead of a 
trial and error procedure. This lends formality to the calibration procedure, decreases 
the arbitrariness and gives measure of likely error in the estimated value. Secondly, the 
best solution is expected when an abundance of data exists. Fortunately, the unsteady 
flow situation can yield an abundance of data, since simple stage measurements over a 
period of time can give many data points. 

In the present study, a formal approach is developed for estimating the Man- 
nings “n” from the measured transient flow data. The developed methodology is based 
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on a nonlinear optimization technique in which the flow equations are embedded as 
equality constraints. 

1.2 Literature review 


The literature relating to identification and estimation of parameters for un- 
steady open channel flow are sparse. The determination of parameters in a set of partial 
differential equations subjected to specified boundary conditions and some measure- 
ments is generally known as the “Inverse Problem” . Significant amount of work relating 
to parameter identification in the groundwater has been reported. Solution techniques 
for aquifer parameter identification problem have included, 

1. Quasilinearation used by Yeh and Tauxe (1971), 

2. Linear programming used by Kleinecke (1971), 

3. Multiple-objective linear programming used by Neumann (1973), 

4. Marquaardt’s nonlinear estimation algorithm used by Garay, Haines and Das 
(1976). 

Becker and Yeh (1972) proposed a methodology for identification of parameters 
in unsteady open channel flow using the influence coefficient approach. The parameters 
chosen for identification was the channel roughness coeflficient and exponent of the 
hydraulic radius in the empirical friction slope relation. The influence coefficients define 
the effect of unit changes in the hydraulic radius exponent and the channel roughness 
coefficient on the errors between observed and simulated values of flow variables like 
velocity and depth. The methodology consists of solving for the objective of minimizing 
the sum of squares of errors between observations and simulated numerical solution of 
the governing equations based on some estimated values of parameters. The algorithm 
starts by specifying initial estimates for the parameters in the feasible region. The 
governing equation are solved by the exphcit (staggered) finite-difierence scheme and 
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the function defining the error is evaluated. The parameters are modified based on the 
evaluated objective function and the influence coefficient matrix. This procedure is 
repeated with new influence coefiicients computed for every iteration. This procedure 
may perform well for lesser number of parameters but proves to be a tedious model for 
estimating a large number of parameters, as simulation has to be carried out several 
number of times. Also the stability and convergence criteria have to be satisfied as an 
explicit scheme is used while solving the governing equations. 

Becker and Yeh (1973a) have used unsteady state flow conditions for identi- 
fication of multiple reach channel parameters by using influence coefiicient approach. 
The parameters estimated were the bed slope and the Mannings “n” at every reach. 
Steady state flow conditions are assumed to be attained prior to the initiation of an 
input flood hydrograph. Steady state computation assumed uniform flow in each reach, 
modified by backwater conditions in the neighbourhood of and some distance upstream 
of each junction. The relatively large rate of change of flow conditions near the junc- ; 
tion necessitates a finer spacial interval in that neighbourhood than elsewhere, in order 
to provide sufficient accuracy and to avoid numerical oscillations that may otherwise ! 
arise. Once initial conditions are determined, unsteady flow equations are solved us- i 
ing staggered explicit finite-difference discretization for interior grid points and for the i 
boundary points. The friction parameters for each reach of channel are identified using ; 
the influence coefficient approach in a sequential manner. Hence, it becomes a tedious | 
procedure. 

The applicability of linear programming together with influence coefficient ap- i 
proach for unsteady open channel flow parameter identification is presented in Becker I 
and Yeh (1973b). The Simplex linear programming algorithm is used in this model. | 
The Minmax criteron (Minimizing the maximum of absolute value of the errors) is ! 
combined with influence coefficient technique to yield a linear programming formula- | 
tion suitable for identification of the parameters of unsteady open channel flow. The I 
parameters estimated are the Mannings “n”, exponent of hydrauHc radius and lateral I 
inflow. I 
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Fread and Smith (1978) proposed a methodology for estimating the parameter 
“n” that varies with space and with stage or discharge and applied the methodol- 
ogy to an hypothetical ideal channel, a mainstem river and a dendritic river system. 
This method uses a gradient type modified Newton-Raphson algorithm to improve the 
initial trial function, such that difference between observed and computed stages is 
minimized. The computed stages are obtained from a four-point implicit finite dif- 
ference solution of one dimensional unsteady flow equations subject to upstream and 
downstream boundary conditions of observed discharge and.stage hydrograph respec- 
tively. 

Wasantha Lai (1995) used the Singular value decomposition method to cali- 
brate the Mannings roughness coefficient in a one-dimensional unsteady flow model. 
The method is used to solve for the parameters after formulating the calibration prob- 
lem as a generalized linear inverse problem. Calibration was repeated with different 
numbers of parameter groups to determine the variation of the output error and un- 
certainity of the parameters with the parameter dimension. Like the earlier methods 
of parameter identification, this method also uses the influence coefficient approach to j 
improve the assumed values of the parameters which are being estimated. However, it | 
employs the Singular value technique so that underdetermined, overdetermined, even 
or mixed determined problems could be solved. Using this approach, optimal grouping 
of sensitive parameters can be determined. I 

In all these previously mentioned works, the parameters are estimated by solv- | 
ing the governing equations in a separate solver for estimated values of the parameters, I 
outside the optimization model or algorithm. The aim of the present study is to esti- j 
mate the parameter “n” by directly embedding the finite-difference discretization of the 
.system of governing equations for the open channel flow into the nonlinear optimization | 
model. As the unsteady flow equations are a set of nonlinear partial differential equa- 
tions, using a nonlinear optimization technique for solution is justified. In this study, j 
the flow equations are discretized using an implicit finite-difference scheme (Preissmann 
scheme). The advantage of using this scheme is that the stability and convergence cri- I 
teria need not be checked. The nonlinear optimization technique used in this study ! 
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is the Sequential Quadratic Programming (SQP), which is available in Numerical Al- 
gorithm Group Library (NAG, 1990). The main advantage of using an embedded 
optimization model for solving the inverse problem is that an iterative procedure is not 
necessary. The flow hydraulics is directly simulated within the optimization model by 
using the finite-difference equations as equality constraints. In the iterative method, 
when the simulation is outside the optimization model, convergence to an optimum 
solution may prove to be difficult if a large number of parameters are to be identified. 
Also the embedding approach results in an elegent formulation and has the potential 
of incorporating various complex flow systems and boundary conditions. However, the 
computational limitations may hinder its universal application. 

In this study, the flow equations are discretized using an Implicit finite-difference 
scheme (Preissmann scheme). The advantage of using this schemeis that the stability 
and convergence criteria need not be checked. 


1.3 Objective of the study 

Objectives of the present study are: 

• To develop an embedded nonlinear optimization model for estimating the rough- 
ness coefficient (Mannings “n”) of an open channel using transient flow observa- 
tion data. 

• To test the developed model for solvability using simulated data as observed data 
for a hypothetical channel. 

• To evaluate the performance of the developed model through illustrative exam- 
ples. 
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1.4 Organization of the thesis 

The thesis is divided into four chapters: 

1. Chapter I contains the introduction to this study and a discussion of the 
relevent literature. 

2. Chapter II presents the methodology and also discusses the optimization 
algorithm used. 

3. Solution results and performance evaluations are dicussed in Chapter III. 

4. Conclusions are presented in Chapter IV. 



Chapter 2 


METHODOLOGY 


2.1 Introduction 


The proposed model for estimating the Mannings roughness coefficient “n” involves 
the use of a nonlinear optimization model with the governing equations of the flow 
forming the equality constraints. Unsteady shallow water flow equations representing 
the mass and momentum conservation are the governing equations. These are a set 
of nonlinear hyperbolic partial differential equations which are discretized using the 
implicit Preissmann flnite-difference scheme (Chaudhry 1993). This scheme is used in 
this study because large computational time steps can be used in the discretization. 
The optimization algorithm used in this work is the Sequential Quadratic Programming 
(SQP) which is available in the NAG Library (NAG, 1990) as subroutine EOdUCF. It 
is designed to minimize an arbitrary smooth nonlinear objective function subjected to 
constraints, which may include simple bounds on the variables, linear constraints and 
smooth nonlinear constraints separately. Since the objective functions and constraint 
functions are smooth i.e., at least twice differentiable, SQP can be used. 
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2.2 Governing Equations 


One-dimensional, unsteady shallow water flow equations for a prismatic, wide rectan- 
gular channel are used as the governing equations in this study. It is however possible 
to extend the methodology to general one-dimensional flow in irregular channels. The 
governing equations adopted are based on the following assumtions: 


• The pressure distribution is hydrostatic. 

• The channel bottom slope is small, so that the flow depths measured normal to 
the channel bottom and measured vertically are approximately the same. 

• The flow velocity over the entire channel cross-section is uniform. 

• The channel is prismatic. 

• The head losses in unsteady flow may be simulated by using the steady-state 
resistance laws, such as the Mannings equation. 


The governing equations represent the mass and momentum conservation and are writ- 
ten as: (Subramanya 1984; Chaudhry 1993) 

Continuity equation 


^ + ^ = 0 
dt dx 


(2.1) 


Momentum equation 


dq d ,q^ 
dt dx^ y 


= gyiso - Sf) 


where: 

y{x, t)=flow depth, 

q{x, t)=discharge per unit width of channel, 
So=slope of channel, 
s/=friction slope, 

x=distance along the flow direction 


( 2 . 2 ) 
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t=tinie. 

The friction slope, s/ is given by the Mannings equation. 


Sf = 


y » 


(2.3) 


in which, n = Maxinings roughness coefficient. 

In the discretized form, the partial derivatives and other coefficients are approximated 
using the Preissmann scheme (Chaudhry, 1993). 


2.2.1 Preissmann Scheme for Discretization 

Equations (2.1) and (2.2) have to be discretized using a finite-difference grid (Fig. 
2.1) before they are used in either a simulation model or in an optimization model as 
equality constraints. 

In the following discussions, subscript i refers to the grid point in the x-direction and 
the superscript k refers to the time index. The partial derivatives and other coefficients 
are approximated as follows, using the Preissmann scheme (Chaudhry, 1993). 


df .. Ut*' + /&') - (ft + fD 

dt 2At 

( 2 . 4 ) 

a/ _ aiftff - ft*') (1 - a)(/i, - ft) 

dx Ax Ax 

( 2 . 5 ) 

f = + ft *') + (1 - + /f)i 

( 2 . 6 ) 


in which, 

Aa:=grid spacing, 

At=computational time step and 
Q!=weighting coefficient. 

In the partial derivatives (Equations 2.4, 2.5 and 2.6), f refers to both q and y: f as a 
coefficient stands for Sf and q. 

The scheme is stable if 0.5 < n < 1.0. A value of 0.8 is used in the present study. 
Substitution of the prece'ding finite difference approximations into Eqs.(2.1) and (2.2) 
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and subsequent rearrangement leads to the followings algebraic equations valid at any 
node i, with ”N” being the total number of nodes being considered. 

Continuity Equation 




2At 


Ax 


+ 


Ax 


(2.7) 


for i = 1 to N - 1 


Momentum Equation 


2At 


+ 


a 

Ax 


riMl 
^ htt\ 


M3!l 4- (l-Q) r (gm)^ _ Ml!l 4. 

^ Ax ^ 




2 Ax' 


2Ax 


{htiir 


+ 
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(A?*' + hi*') + 


+ h'^)] + 


g{l - a)n^ 


( 2 . 8 ) 


for i = 1 to N - 1 


2.3 Formulation of the Problem 


The Inverse problem is formulated as nonlinear optimization model, with the 
resistance coefficient (n) as an unknown variable. The model consists of the open 
channel flow governing equations, boundary conditions and observed values of flow 
variables (depth and discharge). The governing equations are the continuity and mo- 
mentum equations. These governing equations are included as equality constraints in 
the non-linear optimization algorithm, in the discretized form as shown in eqn.(2.7) and 
eqn.(2.8). The continuity equation at each node forms the linear equality constraints 
and momentum equation forms the non-linear equality constraints. The upstream 
boundary conditions are the specified discharge hydrograph at the upstream end. The 
downstream boundary conditions follow the steady-state resistance laws. Hence the 
upstream boundary conditions form the linear constraints and downstream boundary 
conditions form the nonlinear constraints. The objective function is sum of squares of 
the differences between the estimated and measured values of the flow variables. The 
formulation of the Inverse problem can be summarized as shown below: 

Objective Function: The objective function minimizes the sum of squares of the 
differences between the observed flow vaxiable values and the model simulated flow 
variable values. 

Min 


x{M) ^ 


(2.9) 


Where: 

x(0) are the simulated depth or discharge values 
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x(M) are the measured depth or discharge values 

Linear. Constraints: The finite difference forms of the continuity equation and the 
upstream boundary conditions constitute the linear equality constraints of the model. 
These constraints are: 

{yi+i+yi^")-(.yKi + yi) , , (i-Q=)feym 

2Ai Ax Ax 

(2.10) 

for i = lto iV— 1; k = l top (p = Total number of time steps) 

= Specified by inflow hydrograph-, for k = ltop (2-11) 

Non-Linear Constraints: The finite difference form of the momentum equation 
and the downstream boundary conditions constitute the nonlinear equality constraints 
of the model. These constraints are: 


2At 


+ 


a 

Ax 


AMI 

^ httl 


hf+i ^ Ax ^ 


M)! 

h^ 


] + 


gcx 

2Ax 






2Ax 




VH+l ) 


{ht^)f 

r hi+i(qtlf 


] - gsA^ihiit + /if«) + + 1*?)] + 

+ ] — 0/ori = ItoN — l',k = Itop 

ihl)T 


g{l - a)n^ 
2 


( 2 . 12 ) 


for k = 1 top 

In addition, the following physical constraint is specified: 


(2.13) 


n > 0 (2-14) 

The above optimization model is nonlinear because of the nonlinear objective function 
and nonlinear constraints. 
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2.4 Optimization Algorithm 


The above optimization model is nonlinear because of the nonlinear objective 
function and nonlinear constraints. The nonlinear optimization model is solved using 
the Sequential Quadratic Programming (SQP) algorithm, as coded in the NAG library 
(NAG, 1990). The salient features of this algorithm as described in the NAG Users 
Manual is presented here. 

The SQP algorithm is used to solve the general optimization problem: 


Min : F{x) 


(2.15) 


Subject to :l < 


X 

< Alx 

. C(a:) 


< u 


(2.16) 


where F(x) is the nonlinear objective function, Al is an n£, x n coeflhcient matrix, C(x) 
is an n^v element vector of nonlinear constraint. The objective function and constraint 
functions are assumed to be smooth i.e., atleast twice-continuously differentiable. 


At the solution of eq(2.16), some of the constraints will be active, i.e., satisfied exactly. 
An active simple bound constraint implies that the corresponding variable is fixed at 
its bound, and hence the variables are partitioned into fixed and free variables. A 
point X is a first-order Kuhn-Tucker point for eq(2.16) if the following conditions hold 
(Powell, 1974): 


1. X is feasible; 

2. there exists vectors ^ and A (the Lagrange multiplier vectors for the bound and 
general constraints) such that 

j=Cl'A + f (2.17) 

wh6r6 g is th6 grsidisiit of F ovsilusttBd 2 tt x, and == 0 if tho jth. variablo is frsc. 
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3. The Lagrange multiplier corresponding to an inequality constraint active at its 
lower bound must be non-negative, and non-positive for an inequality constraint 
active at its upper bound. 

There are several ways of organizing the matrbc calculations of active set algorithms 
for quadratic programming. By working with the upper triangular matrix method for 
quadratic programming, it is possible to carry out the matrix operations economically. 
The advantage of this approach (Powell, 1983) is that there is a substantial saving of 
computer storage, because the orthogonal matrix is not required. All these features 
are included in E04UCF subroutine (NAG, 1990), which is used in this study. The 
basic structure of E04UCF involves major and minor iterations. The major iterations 
generate a sequence of iterates Xk that converge to z*, a first-order Kuhn-Tucker Point 
of eq(2.16). At a typical major iterations, the new iterate x is defined by 

x = x-+ap (2.18) 

where x is the current iterate, 'a' the steplength and p is the search direction. Also 
associated with each major iteration are estimates of the Lagrange multipliers and a 
prediction of the active set. 

The search direction p is the solution of a quadratic programming subproblem of the 
form. 


Minimize: g^p + \p^Hj, 


Subject toil < { 


P 

AlP 


> < u 


(2.19) 


I Anp J 

where g is the gradient of function at x, the matrix H is a positive definite quasi-Newton 
approximation to the Hessian of the Lagrangian function. An is the jacobian matrix 
of C evaluated at x. The matrix in eq(2.18) is a positive-definite quasi-Newton ap- 
proximation to the Hessian of the Lagrangian function. (Dennis and Schnabel, 1981). 
At the end of each of major iteration, a new Hessian approximation H is defined as a 
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rank-two modification of H. 

where s = x — x (the change in x) 

In E04UCF, H is required to be positive-rdefinite. If H is positive-definite, H defined 
by eq(2.13) will be positive-definite if and only if y'^s is positive. Ideally, y in eq(2.20) 
would be taken bs y^, the change in gradient of the Lagrangian function 

yL=- 9 — ( 2 . 21 ) 

Where denotes the QP multipliers associated with the nonlinear constraints of the 
original problem. If yj^s is not sufficiently positive, an attempt is made to perform the 
update with vector y of the form 

y = yL + ^t^i(ai(x)ci(5) -ai(x)ci(rr)) (2.22) 

where uji > 0. If no succh vector can be found, the update is performed with a scaled 
yx,. After computing the search direction, each major iteration proceeds by determin- j 
ing a steplength a in eq(2.18) that produces a ’sufficient decrease’ in the augmented 
Lagrangian merit function 

L{x, A, s) = F(a:) - k{ci{x) - Si) + I3Pi(ci(a;) - s^)^ (2.23) ! 

where x, A, and s vary during the linesearch. The summation terms in the above ’ 
equation involve only the nonlinear constraints. The vector A is an estimate of the 
Lagrange multipliers for the nonlinear constraints of eq(2.9). Therefore each iteration 
includes, I 

I 

1. The solution of a quadratic programming subproblem. ; 

I 

2. A line search with an augmented Lagrangian merit function and i 

! 

I 

3. A quasi-Newton update of the approximate Hessian of Lagrangian function | 




Chapter 3 


PERFORMANCE EVALUATION 
OF THE PROPOSED MODEL 


3.1 Simulation of Flow Measurement Data 


Performance of the developed optimization model for the estimation of the 
roughness coefficient was evaluated using illustrative example problems for a hypo- 
thetical open channel. In this study, the parameter “n” is estimated for a single reach 
channel and multiple reach channels. As the main aim is to verify the capability of the 
model to estimate the roughness coefficient, without going into the broader issues of 
measurement errors and modelling uncertainities, which are not incorporated. In the 
illustrative examples, identical initial and boundary conditions are implemented in the 
simulation and optimization models. Flow and discharge measurement data in space 
and time were simulated using a numerical simulation model for unsteady open channel 
flow. The Preissmann. discretization scheme was utilized and the solution scheme as 
given by Chaudhry (1993) was adopted for simulation. Using the simulated observation 
data, simulated with known values of “n” also ensures that the estimated “n” values 
and the actual “n” can be correctly compared. 
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3.2 Illustrative Problem 

A rectangular channel with a constant bed slope of 0.0004 is considered. Initial 
steady-state for a single reach channel are given as: y = 2.56m and q = 4m^ /sec/m 
width of the channel. The reach length is 121an and observation stations are assumed 
to be located at 21an interval from the upstream boundary. An incoming triangular 
hydrograph is assumed at the upstream end and observations of discharges and depths 
are simulated by direct numerical solution of eqs(2.1) and (2.2) using the Preissmann 
Implicit finite-difference scheme. The simulation is carried out for three time steps and 
the results are tabulated at each grid locations. There are six nodes for each of the 
time steps as shown in Fig.(3.1). The grid locations are indicated in the Fig.(3.1). The 
succeeding sections discuss the performance evaluation results. 
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Figure 3.1: Grid identification scheme 
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3.3 Verification of the Model for a Single Reach 
Channel 


The ability of the developed model to correctly estimate Mannings “n” for a 
single reach channel is first established. The flow observation data for unsteady flow 
in a channel with a single representative value of “n” is simulated first. The value of 
“n” is assumed to be 0.023 in the simulation model. These simulated observations are 
used as inputs to the optimization model. 

The discretized form of flow equations [eq(2.7) and eq(2.8)] for each box formed by four 
nodes (Fig. 3.1), results in eighteen linear and eighteen non-linear constraints respec- 
tively. The boundary conditions as described below are also needed to be specified. 
The Model solves these flow equations, for given boundary conditions and the flow 
measurments data (here the simulated data), to estimate the parameter “n” (rough- 
ness coefficient). 

Initial Conditions: 

The initial conditions are y(i) = 2.56m and q(i) = 4.0m^/sec for all i = 1 to 7. These 
initial conditions are directly incorporated into the discretized equations (eqs. 2.7 and 
2.8) for the time level k = 1. 

Upstream Boundary Conditions: 

The discharge hydrograph shown in Fig. 3.2 forms the upstream boundary conditions. 
This translates into following linear equality constraints. 


. a:(19, 2) = 4.6667m^/ sec. 
x{2Q, 2) = 5.3333mVsec. 
x(21, 2) = 6.0000m^/sec. 

Downstream Boundary Conditions: 

The downstream boundary conditions are the steady-state resistance laws. 


a:(6,2) 


a:(6, 1)3S02 


(3.1) 

(3.2) 

(3.3) 


n 


(3.4) 
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Figure 3.2: Incoming flow hydrograph at the upstream end of the channel 
(for single reach channel) 


a:(12,2) 
x(18, 2) 


2;(12, 1)3so5 
n 

x(18, 1)3S02 
n 


(3.5) 

(3.6) 


Where: 

x(19,2), x(20,2), x(21,2) = Represents the discharge at grids 19, 20 and 21. 
x(6,2), x(12,2), x(18,2) = Represents the discharge at grids 6,12 and 18 
x(6,l), x(12,l), x(18,l) = Represents the depth at grids 6, 12 and 18 
so = Bottom slope of the channel 
n = Mannings roughness coefficient 

The Objective function is the sum of squares of the difference between the model 
simulated values of the flow variables for the optimal estimate of “n” and the measured 
data of depth and discharge at different nodes. The depth and discharge are treated 
as variables along with the parameters to be estimated i.e., the Mannings “n”. Some 
initial estimate of all these variables are specified to start the optimization algorithm. 
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Embedding technique may prove to be an elegent method for Inverse problems, 
especially when there is abundance of data. Fortunately unsteady flow conditions yield 
numerous data points. To study the robustness of the model, different sets of measured 
data are used. Table. (3.1) presents the simulated measurement data for a single reach 
channel. Actual value of Mannings “n” used in the measurement simulation model is 
0.023. The six scenarios chosen are presented in Table.(3.2). From the results of the 
estimated parameter, it is clear that as number of data points are less, the estimated 
parameter is questionable. Also the optimum value of discharge and depth are almost 
identical to that of the simulated measurement data, when large number of data points 
are input in the optimization model. [Refer to fig.(3.3) and fig.(3-4)]. For scenarios 1 to 
3, the estimated and actual values of “n” are identical. However, when much smaller 
number of measurement data points are used, scenarios 4, 5 and 6, the maximum 
error in the estimation of “n” is about 15 percent only. These results estabflsh the 
correctness of the model solution at least for the simple and ideal case of a single reach 
channel with no measurement errors. 


Table 3.1: Simulated observation Data for a Single Reach Channel 
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Table 3.2: Scenarios Studied for Single Reach Channel 


Scenario No. 

No. of Data Points 

for Depth 

No. of Data Points 

for Discharge 

Estimated 
value of 
“n” 

Actual 
Value of 
“n” 

1 

21 

21 

0.023 

0.023 

2 

17 

17 

0.023 

0.023 

3 

12 

12 

0.023 

0.023 

4 

5 

5 

0.02313 


5 

21 

0 

0.0262 


6 

0 

21 

0.0232 




Figure 3.3: Comparison of flow depths as observed (Simulated) and obtained 
as solution to the optimization model (Single reach, Scenario:!) 

3.4 Estimation of Mamnings “n” for a Multiple Reachl 
Channel 


The formulation is similar to the earlier case of a single reach. The difference 
is that the Mannings “n” to be estimated at each reach is incorporated as a \-ariable 
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Figure 3.4; Comparison of discharge as observed (Simulated) and obtained 
as solution to the optimization model (Single reach, Scenario:!) 

in the optimization model. The model is tested for different sets of simulated flow 
measurement data (depth and discharge). A rectangular channel with six reaches, 
each having the same slope of 0.0004 and with different Mannings “n” is chosen as an 
illustrative example for estimation of “n” values for multiple reaches. The Mannings 
“n” used for simulation at the six reaches of the channel are 0.024, 0.025, 0.021, 0.018, 
0.019 and 0.023 respectively. The initial conditions used in this case are y(i) = 5.56m 
and q(i) = 20.0m^/sec. The discharge hydrograph is shown in fig.(3.5). Table. (3.4) 
presents the simulated data for a Multiple reach Channel. The model is again tested for 
different measurement data and scenarios as shown in tables 3.5 and 3.6 respectively. 
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Table 3.3: Notation used in fig. (3.3) and fig. (3.4) 


Notation used 

Represents for a Single Reach Channel 

esdl 

Depth estimated for 1st time step 

esd2 

Depth estimated for 2nd time step 

esd3 

Depth estimated for 3rd time step 

ssdl 

Simulated measured depth for 1st time step 

ssd2 

Simulated measured depth for 2nd time step 

ssd3 

Simulated measured depth for 3rd time step 

esdil 

Discharge estimated for 1st time step 

esdi2 

Discharge' estimated for 2st time step 

esdi3 

Discharge estimated for 3st time step 

ssdil 

Simulated measured discharge for 1st time step 

ssdi2 

Simulated measured discharge for 2nd time step 

ssdi3 

Simulated measured discharge for 3rd time step 


Table 3.4: Simulated observation data for a Multiple reach Channel 


Grid Point 

Depth 
in metres 

Discharge 
/sec/m width 

Grid Point 

Depth 
in meters 

Discharge 
m^/sec/m width 

1 


20.94 

2 

5.96 

20.52 

3 

5.92 

20.25 

4 

6.09 

19.81 

5 

5.92 

19.33 

6 

6.33 

18.81 

7 

6.74 

23.32 

8 

6.45 

23.14 

9 

6.45 

22:85 

10 

6.59 

22.61 

11 

6.57 

22.33 

12 

6.95 

22.01 

13 


24.79 

14 

6.75 

24.63 

15 

6.77 

24.48 

16 

6.95 

24.30 

17 

6.92 

24.13 

18 

7.31 

23.96 

19 

6.79 

21.67 

20 

7.02 

23.33 

21 

7.33 

25.00 





3.5 Discussion of results 


The evaluation results for multiple reach channels show that the estimated 
values of “n” are fairly accurate for scenarios 7 and 8, where large number of measure- 
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Figure 3.5: Incoming flow hydrograph at the upstream end of the channel 
(for multiple reach channel) 

ment data points are used. For scenario 12, when only 16 data points are used, the 
estimation is not satisfactory. Also, bounds specified on “n” values, are artificially large 
enough to eliminate the influence on the solutions. These evaluation results indicate 
that accurate parameter estimation can be made, if the number of data representing 
measured flow variables are available for a number of spatial locations and time peri- 
ods. This requirement is relevant especially when many parameters are to be estimated 
simultaneously. Since this model uses the embedding approach, it has the potential 
applicability to find other parameters such as bottom-slope, exponent of hydraulic ra- 
dius etc., by treating them as variables in the optimization model. The initial estimate 
of “n” input to the optimization model was found to have insignificant effect on the 
solution results, although the global optimality cannot be guaranieed. 
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Table 3.5: Scenarios Studied for Multiple reach Channel 


Scenario No. 

No. of data points 
of depth 

(y)- 

7 

21 

8 

17 

9 

10 


No. of data points 
of discharge 

( q ) 

Estimated value 
of 

“n” at 
different 
reaches of 
channel 

Actual value 
of 

“n” at 
difierent 
reaches of 
channel 





0.02486 

0.024 


0.02245 

0.025 

21 

0.01887 

0.021 


0.01808 

0.018 


0.01955 

0.019 


0.02195 

0.023 


0.02338 

0.024 


0.02368 

0.025 

17 

0.0191 

0.021 


0.01808 

0.018 


0.0201 

0.019 


0.02185 

0.023 


0.02359 

0.024 ■ 


0.02301 

0.025 

10 

0.01996 

0.021 , 


*0.01893 

0.018 


0.02095 

0.019 r 


0.02159 

0.023 . 
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Table 3.6: Scenarios Studied for Multiple reach channel (contd.) 


Scenario No. 

No. of data points 
of depth 

(y) 

No. of data points 
of discharge 

(q) 

Estimated value 
of 

“n” at 
different 
reaches of 
channel 

Actual value 
of 

“n” at 
different 
reaches of 
channel 




0.02376 

0.024 




0.02219 

0.025 


21 

0 

0.01866 

0.021 




0.01757 

0.018 




0.01857 

0.019 




0.02185 

0.023 




0.02190 

0.024 




0.02579 

0.025 

11 

0 

21 

0.01711 

0.021 




0.02089 

0.018 




0.01671 

0.019 




0.02174 

0.023 




6.89E-06 

0.024 1 




3.76E-02 

0.025 1 

12 

8 

8 

0.01828 

0.021 




0.01770 

0.018 




0.02005 

0.019 




0.02254 

0.023 







Depth In m 
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Figure 3.6: Comparison of flow depth as observed (Simulated) and obtained 
as solution to the optimization model (Multiple reach, Scenario:?) 




Discharge in m3/sec/m width 
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Distance in km 


Figure 3.7: Comparison of discharge as observed (Simulated) and obtained 
as solution to the optimization model (Multiple reach, Scenario:7) 
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Table 3.7: Notation used in fig. (3.5) and fig.(3.6) 


Notation used 

Represents for a Multiple Reach Channel 

emdl 

Depth estimated for 1st time step 

emd2 

Depth estimated for 2nd time step 

emd3 

Depth estimated for 3rd time step 

smdl 

Simulated measured depth for 1st time step 

sind2 

Simulated measured depth for 2nd time step 

smd3 

Simulated measured depth for 3rd time step 

emdil 

Discharge estimated for 1st time step 

eindi2 

Discharge estimated for 2st time step 

emdiS 

Discharge estimated for 3st time step 

smdil 

Simulated measured discharge for 1st time step 

smdi2 

Simulated measured discharge for 2nd time step 

smdi3 

Simulated measured discharge for 3rd time step 






Chapter 4 


SUMMARY AND CONCLUSION 


A mathematical model incorporating the open-channel unsteady flow equa- 
tions within optimization model is developed for estimating the Mannings “n” in a 
single and multiple reach channel system. The discretization of the flow equations is 
carried out using the priessmann implicit finite-difference scheme. This model has the 
potential to estimate parameters like the bed-slope, exponent of hydraulic radius or 
any other parameter, by including the parameter to be estimated as the variable in the 
embedded optimization model. The performance of the developed model is evaluated 
for a number of scenarios. Some of the important conclusion are: 

1. The potential applicability of the developed model is demonstrated using illus- 
trative unsteady open channel flow problem. 

2. As the number of the data points increases, the estimated value of the parameter i 

“n” is very close to the actual value. | 

3. The embedded model is useful especially to estimate a number of parameters; 
simultaneously, which otherwise may become a tedious task when using the sim-j 
ulation model outside the optimization model. 
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